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SUMMARY  t 

The  accurate  prediction  of  satellite  decay  dates  some  months  or  years 
ahead  remains  one  of  the  most  difficult  and  intractable  problems  of  orbital 
mechanics,  chiefly  because  the  lifetime  depends  strongly  on  the  future  variations 
in  air  density,  which  are  at  present  not  accurately  predictable. 

In  this  paper  simple  graphical  methods  are  presented  for  estimating  the 
future  lifetime  of  an  Earth  satellite  from  its  current  rate  of  decay,  using 
theory  adapted  to  an  atmosphere  with  a realistic  variation  of  density  with 
height.  The  effects  of  the  departure  of  the  Earth  and  atmosphere  from  spherical 
symmetry,  and  the  variations  of  density  with  time,  are  approximated  by  specifying 
correction  factors.  Orbits  which  experience  serious  lunisolar  perturbations  call 
for  numerical-integration  methods,  and  the  uses  of  the  computer  programs  PROD  and 
PTDEC  are  described. 


Despite  the  many  uncertainties,  the  remaining  lifetimes  of  most  satellites 
should  be  predictable  with  an  accuracy  of  ±10%  by  these  methods,  which  are  based 
on  many  years'  experience  in  making  the  monthly  decay  predictions  for  the  RAE 
Table  of  satellites. 
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1 INTRODUCTION 

Most  satellites  decay  because  of  the  unceasing  attrition  of  their  orbits 
by  air  drag,  and  their  remaining  lifetimes  depend  on  the  future  variations  in 
air  density,  which  cannot  yet  be  predicted  accurately.  So  lifetime  prediction 
is  inevitably  an  inexact  art:  the  accuracy  is  rarely  better  than  ±10%,  and  some- 
times considerably  worse.  Although  lifetime  prediction  may  seem  to  be  a losing 
battle,  this  is  no  excuse  either  for  an  ignominious  retreat  or  for  the  use  of 
inaccurate  theory,  which  would  add  significant  bias  errors  to  the  inevitable 
random  errors.  The  present  Report  describes  the  simple  graphical  methods 
developed  at  RAE,  based  on  theory  accurate  to  about  2%,  and  the  numerical- 
integration  methods  used  when  lunisolar  perturbations  are  large  enough  to 
invalidate  the  simple  techniques. 

"Why  make  lifetime  estimates?",  you  may  ask.  The  answers  are  various. 

For  scientific  satellites  the  scheduling  of  experiments  may  need  to  be  changed 
if,  for  example,  the  orbital  lifetime  proves  to  be  six  months  instead  of  the 
twelve  months  planned;  or  the  experimenters  may  have  to  estimate  the  latest 
possible  date  for  firing  an  on-board  rocket  to  avoid  premature  decay,  or  decide 
whether  a satellite  will  still  be  in  orbit  at  the  time  of  some  intensive  obser- 
vation campaign,  such  as  the  International  Magnetospheric  Study  of  1978  - if 
not,  a replacement  satellite  might  have  to  be  launched.  The  predicted  lifetime 
also  strongly  influences  the  choice  of  satellites  for  observation  and  orbital 
analysis  over  their  complete  life:  the  amount  of  analysis  and  observing  required 
is  very  different  if  the  lifetime  is  five  years  rather  than  two  years.  When 
a satellite  is  approaching  decay,  the  date  of  decay  needs  to  be  known  in  advance 
if  plans  are  to  be  made  for  intensive  observation  during  the  last  few  weeks  of 
the  life,  for  either  scientific  or  military  motives.  The  need  for  accurate 
decay  predictions  is  particularly  acute  for  optical  observations,  because  the 
periods  of  visibility  for  optical  observation  may  last  for  as  little  as  a week, 
and  a decay  prediction  accurate  to  a few  days  is  needed  to  decide  whether  the 
satellite  will  be  visible  or  unobservable  on  the  day  of  decay. 

Satellites  which  execute  orbital  manoeuvres  on  command  are  of  course  out- 
side the  scope  of  this  paper. 

Lifetime  prediction  is  a thankless  task:  if  your  prediction  is  right, 
others  wrongly  think  it  is  easy;  and  if  you  are  wrong,  you  earn  even  less 
credit.  Our  commitment  to  the  subject  began  in  1958  with  the  compilation  of 
the  first  RAE  Table  of  satellites  recording  details  of  each  successful  new 
launch,  including  estimates  of  orbital  lifetimes.  By  1961  the  Table  had  grown 
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into  a ten-page  RAE  Technical  Memorandum*,  and  from  then  onwards  it  was  issued 

more  frequently:  there  have  been  monthly  issues  for  the  past  twelve  years.  By 

the  end  of  1976  there  were  three  volumes  totalling  more  than  500  pages,  covering 

. 2 3 

about  3000  satellites  and  rockets  in  detail  ’ . The  lifetime  estimates  for 
these  satellites  and  rockets  as  they  near  decay  have  been  updated  at  monthly 
intervals  by  Doreen  Walker  and  the  author,  using  methods  which  have  to  be  very 
simple,  because  a great  many  estimates  must  be  made  quickly  each  month.  In 
developing  methods  for  decay  prediction,  the  main  interest  is  in  lifetimes  of 
a few  months  or  a few  years:  lifetime  estimates  of  more  than  twenty  years  are 
critically  dependent  on  future  solar  activity,  which  is  at  present  unpredictable. 


In  an  atmosphere  which  does  not  vary  with  time,  the  remaining  lifetime  L 
of  a satellite  in  an  orbit  of  eccentricity  e may  be  expressed  in  terms  of  the 
mean  motion  n and  its  rate  of  change  n . The  simplest  theory  gives 


L 


3en 

4n 


(1) 


for  eccentricities  between  0.03  and  0.2.  This  was  the  equation  first  used  in 

4 

1957,  but  improved  versions  were  soon  derived  , and  also  equations  applicable 
for  lower  and  higher  eccentricities,  and  for  an  oblate  atmosphere^  '. 

Strangely  enough,  despite  the  huge  advances  in  space  science  and  techno- 
logy in  the  past  twenty  years,  it  is  still  best  to  use  an  equation  of  this  type, 
with  an  observational  value  of  n , because  the  masses  and  sizes  of  about  90% 
of  newly-launched  objects  are  not  revealed  by  the  launching  authorities.  In 
1976  this  information  was  not  available  on  the  99  Russian  launches  and  the  12 
US  military  launches;  most  of  the  remaining  17  launches  were  of  satellites  which 
entered  high  orbits,  and  lifetime  estimates  were  not  of  much  interest. 

Although  most  satellites  decay  because  their  total  energy  is  whittled  away 
by  air  drag  over  months  or  years,  there  is  an  important  class  of  orbits  - those 
with  eccentricities  between  0.7  and  0.8  - for  which  the  lifetime  is  controlled 
by  the  perturbations  in  the  perigee  height  due  to  lunisolar  gravitational 
attraction.  These  perturbations  can  alter  the  perigee  height  by  more  than 
1000  km  in  a few  years,  and  can  bring  the  perigee  down  to  heights  near  100  km, 
where  rapid  or  catastrophic  decay  occurs.  The  most  important  satellites  in  this 
class  are  the  Russian  Molniyas  and  their  rockets,  of  which  more  than  100  have 
been  lauched.  For  these  satellites  we  rely  upon  numerical  integration  of  the 
lunisolar  perturbations,  backed  up  by  approximate  theoretical  methods.  The 
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numerical-integration  methods  could  of  course  be  used  for  all  satellites,  but 
the  computer  time  required  would  be  prohibitive. 

Sections  2 to  4 of  this  Report  deal  with  the  estimation  of  lifetime  by 
quick  and  simple  graphical  methods,  for  'normal'  decay  under  air  drag.  Section  5 
deals  with  high-eccentricity  orbits  subject  to  strong  lunisolar  perturbations, 
and  describes  the  use  of  the  computer  programs  PROD  (drag-free)  and  PTDEC  (with 
drag  included) . 

2 DEVELOPMENT  OF  APPROXIMATE  THEORY 
2. I Assumptions 

In  developing  theory  for  the  effect  of  air  drag  on  satellite  orbits,  it  is 
simplest  to  assume  a spherical  Earth  with  a spherically  symmetrical  atmosphere 
which  does  not  vary  with  time,  and  has  a density  p varying  exponentially  with 
height  y , so  that 

p = pQ  exp|-(y  - yQ)/H}  , (2) 


where  suffix  0 denotes  a convenient  reference  height  (usually  taken  as  perigee 
height)  and  H , the  'density  scale  height',  is  taken  as  constant.  If  we  write 
z = ae/H  , where  a is  the  semi  major  axis  of  the  orbit,  the  remaining  lifetime 

of  the  satellite  can,  under  these  simple  assumptions,  be  expressed  in  terms  of 
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n as 


for  0.03  < e <0.2.  Within  these  limits  for  e , the  value  of  z is  usually 
between  3 and  40,  so  that  the  neglected  terms  in  (3)  never  exceed  1%.  Similar 
formulae  can  be  derived’’’^  for  e < 0.03  and  > 0.2. 

In  practice,  equation  (3)  needs  correction  because  (a)  the  density  does 
not  vary  exponentially  with  height;  (b)  the  upper  atmosphere  and  the  Earth  are 
appreciably  oblate;  (c)  the  perigee  height  varies  as  a result  of  perturbations 
from  odd  zonal  harmonics  in  the  geopotential;  and  (d)  the  density  does  not 
remain  constant  but  varies  on  daily,  monthly,  semi-annual  and  10  yearly  time- 
scales,  due  to  solar  activity  and  other  effects  (for  a review,  see  Ref  8). 

These  variations  in  density,  often  by  a factor  of  2 and  sometimes  by  a factor  of 
up  to  10,  make  accurate  estimates  of  lifetimes  impossible.  Nevertheless  it  is 
11  reasonable  to  aim  for  an  accuracy  of  about  10Z  for  lifetimes  of  less  than  one  year. 


a mm  i ■ g 
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The  methods  of  lifetime  estimation  presented  here  make  full  allowance  for 
correction  (a),  while  separate  correction  factors  are  specified  in  section  4 for 
corrections  (b),  (c)  and  (d) . 


V p dv  / 


and  assume  that  H varies 


2.2  Variation  of  H with  height 

If  we  define  the  scale  height  H as 
linearly  with  height,  having  a gradient  dH/dy  = y , it  is  found  from  a recent 
revision  of  Ref  9 that  the  extra  terms  to  be  added  inside  the  curly  brackets  in 
equation  (3)  are 


- "(i  - ii) 


+ 0 (v*  ) 


(4) 


In  the  real  atmosphere  y is  not  constant,  but  a constant  value  appropriate  to 
the  initial  perigee  height  will  serve,  because  the  perigee  height  decreases  by 
only  about  $H  during  the  first  90%  of  the  life  . 

Fig  1 shows  the  variation  of  H and  dH/dy  as  given ^ by  the  COSPAR 
International  Reference  Atmosphere  1972  for  high  and  low  levels  of  solar  activ- 
ity: the  curves  are  for  exospheric  temperatures  of  1100  K and  800  K respect- 
ively, corresponding  to  solar  10.7  cm  radiation  energy  F n 7 of  about  150  and 

-22  -2  -1  U,/ 

80  x 10  Wm  Hz  respectively,  at  the  mean  of  the  day-to-night  and 

semi-annual  variations.  The  'high  solar  activity*  is  approximately  the  level 

during  the  solar  maximum  of  1968-70,  and  the  'low'  is  near  the  level  during 

the  prolonged  solar  minimum  of  1974-76.  Although  Fig  1 runs  from  100  km  up  to 

600  km  height,  the  satellites  for  which  lifetime  estimates  are  required  usually 

have  perigee  heights  between  200  and  400  km:  so,  if  a standard  value  of  y is 

needed,  y * 0. 1 is  a suitable  choice. 

2. 3 Formulae  for  lifetime 

On  adding  the  terms  (4),  we  may  rewrite  equation  (3)  as 


for  e < 0.2  and  z > 3 (e  greater  than  about  0.03).  The  maximum  values  of  the  0 
terms  in  (5)  are  0.012  (when  e = 0.2),  0.014  (when  z = 3),  0.005  (if  H = 70  km), 
and  0.01  (if  y = 0.1):  these  maximum  errors  do  not  occur  simultaneously,  and 
the  total  error  due  to  approximations  in  the  theory  may  be  assessed  as  2%. 
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For  z < 3,  the  appropriate  formula  for  lifetime  is' 


3enI0(z)  i 2 , 

{l  - pJ  + 0(e\  p)| 


4nl 


(6) 


where  J , given  in  Fig  14  of  Ref  9,  decreases  almost  linearly  from  1.0  at  z * 0 
to  0.2  at  z = 1.2  and  is  of  order  0.1  for  1.5  < z < 3.  The  functions  I„  and 


Ij  in  (6)  are  Bessel  functions  of  the  first  kind  and  imaginary  argument,  of 
degree  0 and  1 , and  argument  z 
reduces  to 


For  circular  orbits  (z  -*■  0),  equation  (6) 


L = 


3Hn 

2ah 


(1  ~ p) 


(7) 


For  e > 0.2,  the  lifetime  may  be  taken  as 


L = ■njr'  F (e)  ( 1 - *„)  , (8) 

where  F(e)  is  plotted  in  Fig  4.10  of  Ref  11,  and  the  p term  is  carried  over 
from  equation  (5),  the  p/2z  term  being  negligible  for  e > 0.2  (z  > 40).  The 
theory  for  e > 0.2  in  an  atmosphere  in  which  H varies  with  height  has  not  yet 
been  fully  developed,  but  the  assumption  that  the  correction  factor  is  (1  - Jp) 
is  unlikely  to  cause  significant  errors. 

3 RESULTS 

3. 1 Over  the  full  range  of  e 

In  presenting  the  results  it  is  convenient  to  write 

L = Q/n  (9) 

and  to  give  values  of  Q . Fig  2 shows  how  Q varies  with  e for  perigee 
heights  of  200  km  (full  lines)  and  300  km  (broken  line)  for  high  and  low  solar 
activity  over  the  whole  range  of  e from  0 to  0.8.  Over  most  of  this  range  the 
curves  for  perigee  heights  of  200  and  300  km  very  nearly  coincide,  and  the 
curves  may  be  taken  as  applying  over  the  range  200-300  km.  In  preparing  Fig  2, 
the  values  of  Q are  obtained  from  equations  (5),  (6)  and  (8),  using  values  of 
H and  p from  Fig  1. 

The  lifetime  of  any  satellite  perturbed  primarily  by  air  drag  may  be 

estimated  approximately  by  reading  off  Q from  Fig  2 and  dividing  by  the 

2 

observational  value  of  n (rev/day  ).  It  should  be  remembered,  however,  that 
if  e > 0.4,  the  orbit  may  be  strongly  perturbed  by  lunisolar  gravitational 
attraction:  there  is  no  easy  way  of  guessing  the  magnitude  of  the  lunisolar 
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perturbations,  but  it  is  found  in  practice  that  Fig  2 may  be  useful  up  to 
e = 0.8  for  low-inclination  orbits,  while  for  many  high-inclination  orbits  the 

graph  loses  its  accuracy  for  e > 0.6.  In  Fig  2 the  scale  is  terminated  at 

e = 0.8,  because  lunisolar  perturbations  seriously  affect  nearly  all  orbits  with 
eccentricities  greater  than  0.8. 

The  lower  end  of  Fig  2 cannot  be  read  accurately  enough  for  practical  use, 
and  sections  3.2  to  3.4  describe  more  accurate  and  more  detailed  graphs 
(Figs  3 to  10)  for  the  two  levels  of  solar  activity  and  perigee  heights  between 
150  and  600  km. 

3.  2 For  0 < e < 0.03 

Fig  3 shows  the  values  of  Q for  e < 0.03  and  a range  of  values  of  n , 

when  p = 0. 1 and  H has  the  value  given  in  Fig  1.  When  e < 0.01,  the  value 

of  Q is  more  sensitive  to  perigee  height  (or  n)  than  to  e , and  Q is  plotted 
against  n in  Fig  4 for  fixed  values  of  e between  0 and  0.01,  again  for 
p =0.1.  (It  is  more  useful  to  plot  Q against  n rather  than  against  perigee 
height,  because  the  value  of  n is  always  accurately  known,  whereas  the  perigee 
height  is  often  not  available  - it  is  not  given  explicitly  in  the  USAF  2-line 
elements,  which  are  widely  used  in  calculating  lifetimes.)  In  Figs  3 and  4, 
where  p has  its  standard  value  of  0.1,  Q is  given  by 


3enln(z) 

« - - °-1J> 


from  equation  (6):  z = ae/H  is  calculated  using  values  of  H frcm  CIRA  1972  as 
given  in  Fig  1,  and  the  difference  between  the  curves  for  high  and  low  solar 
activity  is  caused  by  the  differences  in  the  values  of  H . Values  of  J are 
read  from  Fig  14  of  Ref  9. 

The  standard  value  of  0.1  for  p in  Figs  3 and  4,  although  convenient 
when  applying  further  theoretical  adjustments,  is  rather  illogical,  because  the 
values  of  H used  are  realistic,  and  vary  with  height  and  solar  activity. 

Figs  3 and  4 have  therefore  been  re-calculated  and  re-drawn  using  the  values  of 
both  H and  p given  in  Fig  1.  The  results  are  shown  in  Figs  5 and  6,  which 
may  be  regarded  as  giving  values  of  Q based  on  realistic  values  of  both  H 
and  p in  the  best  existing  reference  atmosphere. 

Figs  5 and  6 therefore  provide  the  best  set  of  values  for  Q for  use  in 
practical  estimation  of  lifetime,  though  still  subject  to  the  errors  to  be  dis- 
cussed in  section  4.  The  values  of  Q are  given  in  rev/day  on  the  assumption 
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that  n is  calculated  by  differencing  two  values  of  n (rev/day)  at  least  a 

• 2 

day  apart,  and  dividing  by  the  time  interval,  to  give  n in  rev/day  . The 
lifetime  Q/n  is  then  in  days.  It  is  intended  that  Fig  5 should  be  used  if 
0.01  < e < 0.03,  and  Fig  6 if  0 < e < 0.01. 

Fig  5 shows  that  the  influence  of  solar  activity  on  Q is  greatest  when 
e is  small:  Q is  about  25%  greater  for  high  solar  activity,  for  n < 16.  This 

merely  reflects  the  fact  that  H is  about  25%  greater  when  solar  activity  is 

high,  for  y^  > 250  km  (see  Fig  1),  since  Q is  proportional  to  H for  e = 0 
and  p is  small  for  y^  > 250  km.  Although  Q is  greater  when  solar  activity 
is  high,  the  density  is  much  greater,  and  so  therefore  is  n . Thus  the  life- 
time Q/n  is  always  shorter  when  solar  activity  is  high. 

The  differences  and  peculiarities  in  the  shapes  of  the  curves  in  Figs  3 to 
6 stem  from  the  fact  that  they  embody  the  variations  of  H in  Fig  1,  while  Figs 
5 and  6 are  also  affected  by  the  variations  of  p in  Fig  1.  For  example,  the 
strong  decrease  of  Q in  Fig  4 as  perigee  height  drops  from  300  km  to  150  km  is 

a direct  result  of  the  strong  decrease  in  H , from  52  km  at  300  km  height  to 

18  km  at  150  km  height  (for  high  solar  activity).  This  effect  is  even  greater 
in  Fig  6,  because  of  the  increase  in  p as  perigee  height  drops  from  300  to 
150  km. 

3. 3 For  0.03  < e 4 0.2 

When  e is  within  this  range  of  values,  Q increases  almost  linearly 
with  e and  the  most  accurate  graphical  presentation  is  achieved  by  writing 


and  plotting  Q/e  against  e . The  values  of  Q/e  are  given  by  the  right-hand 
side  of  equation  (5),  with  a factor  3n/4  outside  the  curly  brackets,  rather  than 
3en / (4n) . 

Figs  7 and  8 give  the  values  of  Q/e  for  perigee  heights  between  150  and 
600  km,  for  high  and  low  solar  activity,  with  values  of  H from  Fig  1.  In 
Fig  7 the  value  of  p is  taken  constant  at  0.1,  while  in  Fig  8 realistic  values 
of  p are  used,  taken  from  Fig  1.  Because  of  the  wide  scale,  accurate  values 
of  Q/e  can  easily  be  read  from  Figs  7 and  8.  Since  the  range  of  variation  of 
Q/e  is  so  small,  the  variations  in  p (which  has  a minimum  near  400  km)  are 
enough  to  convert  the  curves  of  Fig  7 to  an  arched  form  (usually  with  maxima 
between  300  and  400  km)  in  Fig  8.  In  practice,  satellites  with  lifetimes  of  a 
few  weeks  or  months  usually  have  perigee  heights  between  150  and  300  km,  and 
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Fig  8 shows  that  the  values  of  Q are  then  confined  to  a very  narrow  band: 
at  e = 0.1,  for  example,  1.10  < Q < 1.15  for  high  solar  activity,  and 
1.11  < Q < 1.15  for  low  solar  activity. 

3.4  For  0.2  4 e ■$  0.8 

The  value  of  Q , now  given  by  equation  (8),  is  plotted  against  e for 
y = 0.  1 in  Fig  9,  for  perigee  heights  of  200,  400  and  600  km:  the  variation  with 
perigee  height  is  linear.  Fig  10  shows  the  values  of  Q when  y has  the  CIRA 
values  of  Fig  1,  for  perigee  heights  between  150  and  600  km.  With  these  real- 
istic values  of  y , there  is  little  variation  of  Q with  perigee  height  for 
perigee  heights  between  150  and  500  km.  When  e > 0.4,  Fig  10  should  be  used 
with  caution,  because  the  orbit  may  suffer  significant  perturbations  in  perigee 
height  due  to  lunisolar  gravitational  effects.  If  these  perturbations  are  less 
than  about  50  km,  it  may  still  be  possible  to  use  Fig  10,  if  the  value  of  n 
can  be  averaged  over  a cycle  of  the  lunisolar  perturbation;  but  if  the  perturba- 
tions are  larger  than  100  km,  numerical-integration  methods  are  needed  (see 
section  5). 


4 CORRECTIONS  TO  THE  RESULTS 
4. 1 Introduction 

Figs  5,  6,  8 and  10  provide  values  of  Q in  an  atmosphere  with  realistic 
values  of  H and  y , for  all  values  of  e and  perigee  heights  between  150  and 
600  km.  These  graphs  would  give  accurate  lifetimes  if  the  Earth  and  atmosphere 
were  spherical  and  the  density  did  not  vary  with  time.  These  ideal  conditions  do 
not  apply  in  practice,  and  several  corrections  are  necessary,  which  will  now  be 
described,  although  they  are  not  easily  presented  in  graphical  form. 


Variation  of  perigee  height 


One  important  correction  arises  because  in  practice  the  perigee  height 
varies  with  the  argument  of  perigee  to  . If  we  write 


(12) 


where  r is  the  perigee  distance  from  the  Earth's  centre  and  R the  local 
P P 

Earth  radius,  the  variation  in  y arises  from  (a)  the  variation  of  r , due 

P P 

mainly  to  the  effects  of  odd  zonal  harmonics  in  the  geopotential,  and  (b)  the 

variation  of  R^  with  latitude.  In  terms  or  the  argument  of  perigee  u , the 

variation  in  r is  of  the  form 
P 

r “ r„n  = ” a6  sin  u + (Ar  ) + (Ar  ) , O- 

p pU  p 2 p LS 


I 
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where  aB  is  the  amplitude  of  the  oscillation  due  to  odd  zonal  harmonics,  given 
in  Fig  11  of  Ref  12,  and  of  order  10  km  for  most  inclinations;  and  suffix  0 

denotes  a reference  value,  chosen  tc  be  where  w = 0.  The  quantity  (Ar  is  the 

- - - - — - p 

amount 


by  which  r differs  from  a(l  - e);  it  is  usually  of  order  1 km  and 
P # 14 

will  be  ignored  here.  (Ar  ) c is  the  lunisolar  perturbation  , which  is  usually 

p Lu 

less  than  1 km  for  orbits  with  e < 0.4  and  will  also  be  ignored:  if  (Ar  ) c is 

p Ltu 

large,  numerical-integration  procedures  have  to  be  used  (see  section  5).  The 

decrease  in  r due  to  air  drag  is  allowed  for  in  the  lifetime  formulae  and 
P 

need  not  be  considered.  The  local  Earth  radius  at  perigee,  R , may  be 

P 

expressed  as 


R 


.2.  .2 

= 6378.1  - 21.4  sin  i sin  u> 


'p  ' yPo 


km  (14) 

Hence,  if  y n is  the  value  of  y when  w = 0,  equations  (12)  to  (14)  give 

pu  p 

2 2 

.n  i sin  w - aB  sin  to  km  . (15) 

The  variation  of  y with  u for  a variety  of  inclinations  i is  shown  in 

P \2 

Fig  11.  (The  values  of  aB  used  in  Fig  11  are  from  Ref  15:  more  recent  values 

are  within  1 km  except  near  i = 10°.) 

The  best  method  of  allowing  for  the  variation  of  yp  with  w is  to  modify 

the  current  value  of  n so  that  it  becomes  the  value,  ^ , say,  that  n would 

have  if  the  current  value  of  y were  equal  to  its  mean  value  y during  the 

P P • 

subsequent  lifetime  (again  excluding  the  decrease  due  to  air  drag) . Since  n 

is  proportional  to  , equation  (2)  gives 


n exp | (y  - y )/h} 


(16) 


P 'P 

where  H is  given  in  Fig  1. 

If  the  argument  of  perigee  is  expected  to  move  through  several  revolutions 
before  decay,  equation  (15),  averaged  over  a cycle  of 


gives 


p0 


2 

+ 10.7  sin  i 


km 


(17) 


so  that,  from  (16), 


n exp | (y  - ypQ  - 10.7  sin2i)/H| 


(18) 

So  the  best  plan  is  to  read  off  (y  - y »)  from  Fig  II  and  H from  Fig  1,  and 

P 

take  the  lifetime  as  Q/fi  , where  n is  given  by  (18).  Alternatively,  to  avoid 
this  correction,  n may  be  evaluated  when  (y^  - y^g)  is  within  1 km  of  its  mean 
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value  10.7  sinl 2 * * *i  (indicated  on  the  right  of  Fig  11).  For  i = 90°,  for  example, 
this  would  mean  that  n should  be  evaluated  when  w is  within  the  ranges: 
68-82°,  98-112°,  208-214°  or  326-332°.  Fig  12  shows  how  these  'safe'  values  of 
u , where  no  perigee  height  correction  is  needed,  vary  with  i : the  'safe  areas' 
are  stippled.  The  safest  values  of  w are  in  the  northern  hemisphere  and 
generally  not  far  from  90°.  For  20°  < i < 55°,  the  best  value  of  w is  90°; 
for  68°  < i < 90°,  a good  simple  rule  is  to  take  w = i - 10°  (or  190°  - i) ; for 
55°  < i < 68°,  the  only  safe  values  of  co  are  two  narrow  bands,  in  the  southern 
hemisphere  near  205°  and  335°  for  55°  < i < 62°. 


If  the  argument  of  perigee  is  expected  to  travel  through  less  than  about 

lj  revolutions  before  decay,  ie  if  Lu  is  less  than  about  540°,  the  mean  value 

of  (y^  - Ypg)  i-n  the  remaining  life  may  be  estimated  graphically  from  Fig  11. 

The  difference  between  the  current  value  of  (y_  - y .)  and  this  mean  value 

P pU 

# . r / p 

(y  - y is  substituted  into  (16)  to  give  n and  hence  L = Q/n  . It  is  not 
P 

always  easy  to  estimate  the  mean  value  of  (y  - y n)  from  Fig  11,  and  it  is  use- 

2 P P^  0 

ful  to  note  that  10.7  sin  i is  the  mean  value  between  u)  = 90°  and  oo  = 270°,  as 

well  as  over  a whole  cycle.  In  practice  an  iterative  procedure  may  be  needed: 

L is  first  calculated  using  (18);  but  if  it  then  turns  out  that  Lu  < 540°,  a 

recalculation  of  y^  is  necessary,  leading  to  a revised  value  of  L from  (16). 

4. 3 Atmospheric  oblateness 


The  correction  specified  in  section  4.2  tacitly  takes  account  of  the  main 
effects  of  atmospheric  oblateness,  because  of  the  underlying  assumption  that  the 
density  depends  on  the  height  y^  above  an  oblate  Earth.  This  is  equivalent  to 
assuming  that  the  constant-density  surfaces  are  oblate,  and  of  ellipticity 
e =*0.0033.  However,  the  theory6  indicates  a further  small  change  in  the  life- 

• l £rn^2^Z^  2 1 

time  formula,  by  a factor  7l  + ^ sin  i cos  2S|  for  z < 5,  and 

l 2e  . 2.  -)  0 

v + Sln  1 cos  2u>  for  z > 5.  Here  cos  2u  is  the  mean  value  of  cos  2u 

in  the  remaining  lifetime.  The  changes  due  i o this  factor  can  be  as  much  as 

10%  for  polar  orbits  with  cos  2u  & 1 and  z =*  5 ; but  these  conditions  are 
most  unlikely  to  occur  simultaneously  and  for  an  'average'  orbit,  with 
i =*  45°  , cos  2u  ^ { and  z =*  10  (or  1),  the  change  is  2%.  In  view  of  the 
much  larger  possible  errors  due  to  variations  in  air  density,  this  correction  is 
rarely  worth  making. 


4.4  The  future  variations  of  air  density 

This  is  by  far  the  largest  source  of  error,  and  a very  intractable  one. 
The  variations  are  considered  here  under  four  headings:  solar-cycle,  day-to- 
night;  semi-annual;  and  irregular  day-to-day. 
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4.4.1  Variations  over  the  1 1-year  solar  cycle 

Fig  13  shows  the  variation  of  density  with  height  by  day  and  by  night  at  a 

g 

typical  sunspot  maximum  and  minimum  . At  heights  near  600  km  the  density  varies 
by  a factor  of  20  during  the  sunspot  cycle,  with  corresponding  changes  in  satel- 
lite decay  rates.  At  present  neither  the  intensity  nor  the  date  of  a sunspot 
maximum  can  be  predicted  accurately,  so  lifetime  estimates  can  only  take  account 
of  solar-cycle  variations  in  a very  approximate  way.  Lifetime  estimates  are 
most  often  made  for  satellites  with  perigee  heights  between  200  and  300  km,  and 
for  these  satellites  the  variation  in  density  during  a solar  cycle  is  less,  the 
ratio  of  maximum/minimum  density  being  between  2 and  5. 

If  the  satellite  is  expected  to  remain  in  orbit  for  several  solar  cycles, 
and  a fairly  accurate  value  of  n is  available,  the  lifetime  L may  be  taken 
as 

(19) 


‘ *(»  ■ 


1 1 1 


where  p is  the  current  density  at  perigee  height,  and  p is  its  mean  value 
P _ P 

over  a solar  cycle.  If  the  variation  of  solar  activity  with  time  during  a solar 

cycle  has  a sinusoidal  form  - admittedly  not  always  a good  approximation  - we 

may  take  p = Up  + p ] , where  'min'  and  'max'  denote  values  averaged 

P ' Fmin  Pmax  ' 

over  day  and  night  at  the  minimum  and  maximum  of  the  cycle.  Values  of  p , 

P 

taken  as  the  arithmetic  mean  of  the  CIRA  values  for  exospheric  temperatures  of 

800  and  1100  K,  are  plotted  as  a broken  line  in  Fig  13:  it  nearly  coincides  with 

the  night-time  solar  maximum  curve.  For  a perigee  height  of  300  km,  the 

correction  factor  p /p  is  between  0.2  and  2;  for  200  km  it  is  between  0.5 

P P 

and  1.5. 

For  a satellite  which  is  not  expected  to  remain  in  orbit  for  more  than  a 
fraction  of  a solar  cycle,  the  density  p^  in  equation  (19)  needs  to  be  the 
mean  density  during  the  remaining  lifetime.  This  is  extremely  difficult  to 
estimate,  particularly  at  solar  minimum,  when  it  is  not  known  how  soon  the 
solar  activity  will  begin  to  increase.  For  example,  in  1975  it  seemed  likely 
that  solar  activity  would  increase  in  1976;  in  1976,  a sharp  increase  early  in 
1977  seemed  almost  inevitable.  But  neither  expectation  was  fulfilled. 

The  only  good  thing  to  be  said  about  the  solar-cycle  variation  is  that  it 
loses  its  importance  for  satellites  with  lifetimes  of  only  a few  months  - and 
these  form  the  most  important  group. 

4.4.2  Satellites  with  lifetimes  of  several  solar  cycles 

If  the  satellite  lifetime  exceeds  about  three  solar  cycles,  the  rate  of 
decay  is  very  slow,  and  it  is  often  difficult  to  calculate  a current  value  of 
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n accurately.  It  is  then  useful  to  express  L in  terms  of  the  raass/area  ratio 
of  the  satellite  m/S  , assuming  a density  profile  averaged  over  a solar  cycle. 
This  average  is  taken  to  be  the  same  as  in  the  past  few  cycles,  the  density 
being  taken  as  that  given  by  the  CIRA  1972  model  with  exospheric  temperature 
900  K:  this  density  is  slightly  higher  than  that  given  by  the  mean  curve  in 
Fig  13.  It  should  be  emphasized,  however,  that  the  average  density  in  the  future 
may  differ  greatly  from  this  value,  particularly  if  in  the  next  50  years  there 
is  a period  of  low  solar  activity  akin  to  the  'Maunder  minimum'  which  occurred 
between  1640  and  1710,  when  sunspots  were  a rarity'^’ 


An  equation  for  L(S/ra)  as  a function  of  perigee  density  p and  eccen- 

P 

tricity  may  be  found  from  equation  (7.6)  of  Ref  11,  namely 


P 


P 


T5(e) 
6 /aH 


(20) 


where  5(e)  is  a function  of  e given  in  Fig  7.1  of  Ref  11.  In  equation  (20) 

• % 2 

the  rate  of  change  of  period  T (days/day)  may  be  written  as  - n/n  ; also 

6 = FSCp/m  , where  the  drag  coefficient  will  for  the  moment  be  taken  as 

2.2,  and  F is  a factor  which  may  be  taken  as  0.91  ± 0.05  for  0 < i < 75° 

(the  error  is  larger  for  i > 75°,  but  still  much  smaller  than  the  possible 
errors  in  the  assumed  future  density  values).  Thus  equation  (20)  becomes 


(21) 


• 2 2 

where  n is  in  rev/day,  n in  rev/day  , m in  kg,  S in  m , H and  a in  m, 

3 — 6 “1  5 

and  p in  kg/m  . Taking  n = 274.5(10  a)  ‘ , and 
_c  P -3  . 

10  a = (6.37  +10  y )/(l  - e)  , where  y is  the  perigee  height  in  km, 

P ' P 

equation  (21)  may  be  rewritten  as 

G(yp)Z(e)  = , (22) 


c / y “ 200 \- 2. 5 

where 

G(yp)  - 10  pp/5(l  * -^m~) 

(23) 

and 

l (e)  - 1.362(1  - «)2,5/3(e) 

(24)* 

* Although  5(e)  depends  primarily  on  e , it  is  also  dependent  on  H when  e 
is  small.  When  e ^ 0.02,  we  may  write  E(e)  as 

12.84(a/H)i  exp(-  z)I0(z){l  - e(2.5  - 21^)}  . 

When  e = 0.02,  the  variation  of  E(e)  is  about  ±4%.  Values  of  E(e)  for 
e > 0.2  are  shown  in  Fig  21. 
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From  equations  (22)  and  (9)  we  have 

L(S/m)  = Q/{365G(yp)i:(e)(  (25) 

if  L is  to  be  in  years  rather  than  days. 

Fig  14  shows  the  values  of  L(S/m)  given  by  equation  (25)  for  perigee 

heights  of  200-2000  km,  and  eccentricities  between  0 and  0.8.  The  values  of 

Q (=  Ln)  are  obtained  from  equations  (5),  (6)  or  (8).  G(y  ) is  obtained  from 

P 

equation  (23)  and  CIRA  1972,  and  E(e)  from  equation  (24)  and  Fig  7.1  of  Ref  II. 
The  density  profile  in  CIRA  1972  leads  to  some  wide  excursions  in  the  values  of 
p at  heights  between  750  and  1000  km,  which  lead  to  corresponding  excursions  in 
the  curves  for  L : these  have  been  smoothed  out  in  drawing  the  curves  in  Fig  14, 
to  assist  in  interpolation.  It  should  be  emphasized  again  (a)  that  the  curves  of 
Fig  14  are  subject  to  considerable  uncertainty,  because  the  future  level  of  solar 
activity  in  the  next  50  to  100  years  is  not  known,  and  (b)  that  the  results  only 
apply  when  the  satellite  experiences  an  average  density  corresponding  to  that 
given  in  CIRA  1972  for  an  exospheric  temperature  of  900  K - the  average  over 
recent  solar  cycles.  This  means  that,  in  general,  the  curves  of  Fig  14  should 
only  be  used  if  the  satellite  is  in  orbit  for  at  least  three  solar  cycles,  or 
approximately  30  years.  However,  the  curves  can  sometimes  give  a good  indica- 
tion of  shorter  lifetimes:  if  a satellite  is  launched  at  solar  minimum,  and 
Fig  14  indicates  a lifetime  of  5 years  or  11  years,  these  values  might  be  quite 
accurate  if  it  turns  out  that  the  5 years  are  from  minimum  to  maximum  of  a 
typical  solar  cycle,  or  the  11  years  make  up  a complete  and  fairly  average 
solar  cycle.  In  these  special  circumstances,  the  density  experienced  by  the 

satellite  will  be  'average'.  A diagram  similar  to  Fig  14  has  been  given  by 

_ ■ • 18 
Penni 

To  obtain  the  lifetime  in  years,  the  values  of  L(S/m)  in  Fig  14  must  be 

2 

multiplied  by  the  mass/area  ratio  of  the  satellite,  m/S  , expressed  in  kg/m  . 

2 

The  value  of  S is  the  average  cross-sectional  area,  equal  to  i "rrd  for  a sphere 

of  diameter  d , and  approximately  equal  to  0.8Zd  for  a tumbling  cylinder  of 

length  l and  diameter  d . The  numerical  value  of  m/S  is  often  near 
2 

100  kg/m  , and  the  values  of  L in  Fig  14  would  then  lie  between  10  years  and 
lO^years.  Values  of  m/S  for  a variety  of  satellites  are  given  in  Fig  15:  the 

O 

values  are  between  50  and  200  kg/m  for  nearly  all  satellites. 

When  the  perigee  height  exceeds  600  km,  the  drag  coefficient  assumed  here, 

19 

2.2,  is  likely  to  increase  , possibly  to  values  as  high  as  3,  because  of  the 
decrease  in  the  ratio  of  satellite  velocity  to  molecular  velocity.  No  allowance 


r 
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has  been  made  here  for  this  increase,  because  its  exact  form  is  not  known;  if 
0^  is  known  to  differ  from  2.2,  the  lifetime  L should  be  multiplied  by  2.2/C^. 

The  results  presented  in  this  section  are  of  course  not  strictly 
'corrections’  to  the  previous  methods,  but  rather  an  alternative  method,  when 
the  usual  n-method  becomes  impracticable,  because  the  current  value  of  n 
cannot  be  accurately  estimated. 

4.4.3  The  day-to-night  variation 

The  day-to-night  variation  in  air  density  is  not  so  large  as  the  solar- 

cycle  variation,  as  Fig  13  shows.  But  the  maximum  daytime  density  can  exceed 

the  minimum  night-time  density  by  a factor  of  up  to  5 at  500  km  height,  so  that 

a lifetime  calculated  from  the  value  of  n when  the  perigee  density  is  at  the 

night-time  minimum  can  be  in  error  by  a factor  of  3 if  the  perigee  later 

experiences  the  average  density.  If,  as  is  more  usual,  the  perigee  height  is 

between  200  and  250  km,  the  error  factor  is  less,  though  still  important  - 

having  a maximum  value  near  1.5.  Usually  perigee  goes  through  a day-to-night 

cycle  in  a few  months.  As  with  the  solar-cycle  variation,  the  remedy  is  to 

adjust  the  lifetime  by  multiplying  by  p / p where  p is  the  current  density 

_ P P P 

at  perigee  and  p^  the  mean  value  over  the  remaining  lifetime,  which  would  be 
the  mean  value  over  a day-to-night  cycle  if  perigee  is  to  sample  several  such 
cycles . 

The  day-to-night  variation  is  somewhat  intractable,  because  even  the 
current  local  time  at  perigee  is  not  readily  available  and  needs  to  be  calcu- 
lated; its  future  variation  is  not  accurately  predictable,  because  it  depends 
on  the  rate  of  decay,  which  itself  depends  on  the  variation  of  local  time. 

For  satellites  of  several  years  life,  the  best  policy  is  to  calculate  n at  a 
time  when  the  density  at  perigee  is  near  the  mean  of  the  day-to-night  variation, 
that  is,  when  the  local  time  at  perigee  is  near  09 h or  20h  local  time. 

The  only  helpful  feature  of  the  day-to-night  variation  is  that  it  becomes 

small  as  decay  approaches:  for  a perigee  height  of  250  km  the  maximum  value  of 

p /p  (for  low  solar  activity)  is  1.4,  decreasing  to  1.2  at  200  km.  When  solar 
P * P 

activity  is  high,  the  corresponding  values  are  1.2  and  1.1.  Also,  when  the 

satellite  is  within  a month  or  two  of  decay,  the  mean  value  of  p during  the 

P 

remaining  life  is  likely  to  be  nearer  to  the  current  value  than  to  the  mean 
value  over  a complete  day-to-night  cycle. 
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4.4.4  The  semi-annual  variation 

In  addition  to  the  wide  variations  shown  in  Fig  13,  the  air  density  under- 
goes a fairly  regular  semi-annual  variation  with  maxima  usually  in  April  and 
October,  and  minima  usually  in  January  and  July.  Although  the  variations  are 
smaller  in  amplitude  than  the  solar-cycle  and  diurnal  effects  - the  factor  of 
variation  is  about  1.6  at  heights  near  200  km,  increasing  to  about  3 at  500  km  - 
the  semi-annual  variation  is  more  important  for  lifetime  estimation  (a)  because 
it  extends  down  to  150  km  height  without  much  diminution,  and  (b)  because  the 
3-month  interval  between  maximum  and  minimum  introduces  considerable  errors 
into  the  most  important  lifetime  estimates,  those  of  between  2 and  6 months. 

The  semi-annual  effect  varies  from  year  to  year  in  amplitude  and  phase,  and 

cannot  yet  be  predicted  accurately.  A standard  variation  is  given  in  CIRA , but 

this  has  proved  to  be  inadequate  in  the  1970s,  and  it  is  preferable  to  use  the 
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variation  of  density  recommended  by  Walker  , shown  in  Fig  16.  These  values,  V , 
are  based  on  results  in  1972-5  at  heights  of  200-250  km,  which  is  the  most 
likely  range  of  perigee  heights  for  satellites  with  lifetimes  of  a few  months. 

To  allow  for  the  semi-annual  variation,  the  lifetime  L should  be  multi- 
plied by  V/V  , where  V is  the  value  at  the  current  date,  and  V is  the  mean 
value  of  V during  the  remaining  lifetime.  Again,  an  iterative  procedure  may 
be  needed  if  V/V  differs  greatly  from  1. 

In  Fig  17  the  semi-annual  correction  V/V  is  plotted  against  lifetime, 
for  lifetimes  of  up  to  160  days,  for  predictions  made  on  particular  days  of  the 
year  - year  days  0,  20,  40,  60  ...  340.  The  correction  for  intermediate  dates 
can  be  found  by  interpolation.  Although  subject  to  some  uncertainty  because  of 
the  yearly  variations  in  the  amplitude  and  phase  of  the  semi-annual  variation. 

Fig  17  provides  an  easy  and  reasonably  accurate  correction,  which  needs  to  be 
used  for  all  lifetime  estimates  of  less  than  160  days. 

Fig  18  shows  the  variation  of  V/V  during  the  course  of  the  year  for 
satellites  with  lifetimes  of  (a)  20-100  days  and  (b)  100  days  to  1 year. 

For  satellites  with  lifetimes  of  many  years,  errors  due  to  the  semi-annual 
variation  may  be  minimized  by  taking  values  of  n calculated  near  the  mean  of 
the  semi-annual  variation,  that  is,  near  19  February,  14  May,  17  September  or 
11  30  December,  where  V/V  « 1 . 
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4.4.5  Irregular  day-to-day  variations 

Air  density  at  heights  above  150  km  suffers  continual  day-to-day  variations 
of  order  ±10%,  and  much  larger  variations  (by  a factor  of  up  to  2 at  200  km  and 
6 at  600  km)  at  times  of  major  geomagnetic  disturbances.  These  irregular  day-to- 
day  effects  must  be  regarded  as  largely  unpredictable  at  present,  although  some 
of  the  likely  variations  may  be  partially  predictable  from  the  27-day  recurrence 
tendency  in  solar  activity:  forecasts  of  solar  activity  and  geomagnetic  distur- 
bances up  to  3 days  ahead  are  made  daily  by  the  USAF  Geophysical  Warning  Center. 
The  skilled  predictor  will  take  advantage  of  any  available  information  of  this 
kind,  which  can  be  most  useful  for  decay  predictions  at  the  very  end  of  a 
satellite's  life.  The  possibility  of  unexpected  magnetic  storms  remains,  how- 
ever, and  predicted  lifetimes  of  about  4 days  may  prove  to  be  too  large  by  a 
factor  of  2 if  a major  geomagnetic  storm  occurs  on  the  day  after  the  prediction 
is  made. 

In  predicting  decays  several  months  ahead,  future  day-to-day  variations 
are  less  important;  but  the  current  rate  of  orbital  decay  may  be  calculated  at  an 
extreme  of  the  current  day-to-day  irregularities.  To  avoid  this  danger,  it  is 
usually  better  to  evaluate  the  current  rate  of  decay  over  at  least  2 or  3 days. 

A subtler  stratagem,  for  an  important  decay,  is  to  make  predictions  every  day  or 
two  and  to  take  a weighted  mean  value,  with  weighting  inversely  proportional  to 
the  remaining  lifetime,  as  proposed  in  Ref  31. 

4.5  Envoi 

After  hearing  of  all  these  corrections,  so  many  of  them  intractable  or 
unpredictable,  anyone  who  contemplates  making  predictions  may  well  feel  like 
giving  up  in  despair  and  may  seek  refuge  in  elaborate  computer  programs  of 
numerical  integration.  This  apparent  escape  route,  discussed  in  section  5,  is 
partly  illusory,  because  the  solar  activity  and  semi-annual  corrections  still 
have  to  be  specified  separately,  and  the  peril  of  the  unpredicted  geomagnetic 
storm  still  lurks  around. 

Despite  the  difficulties  of  making  all  the  corrections,  the  graphical 
methods  do  often  yield  lifetime  estimates  accurate  to  about  10%,  though  much 
larger  errors  can  easily  be  made  if  a necessary  correction  is  neglected. 

5 PROCEDURES  WHEN  LUNISOLAR  PERTURBATIONS  ARE  LARGE 

5. 1 The  PROD  program 

If  the  perigee  height  suffers  large  perturbations  due  to  lunisolar  gravita- 
tional attraction,  the  approximate  analytical  methods  described  in  sections  2 to  4 
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become  ineffective,  and  it  is  necessary  to  resort  to  numerical  integration.  For 

this  purpose  the  computer  program  PROD  (program  for  orbital  development)  was 
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developed  by  G.E.  Cook  in  the  1960s.  To  avoid  excessive  use  of  computer  time, 
the  effects  of  solar  radiation  pressure  and  air  drag  are  omitted  in  PROD  and  the 
lunisolar  disturbing  functions  are  averaged  analytically  with  respect  to  the 
mean  anomaly  of  the  satellite.  Zonal  harmonics  in  the  geopotential  and  harmonics 
in  the  expansion  of  the  solar  and  lunar  potential  up  to  the  ninth  can  be  included, 
although  in  practice  only  the  second  harmonic  of  the  solar  potential  and  up  to 
the  third  in  the  lunar  are  required.  The  interval  of  integration  can  be  chosen 
by  the  user. 

PROD  has  been  applied  in  a large  number  of  long-term  lifetime  predictions. 

It  has  been  particularly  valuable  in  predicting  the  lifetimes  of  Molniya 
satellites  with  eccentricities  near  0.7  initially.  The  usual  effect  of  luni- 
solar perturbations  on  these  satellites,  and  on  many  rockets  left  in  transfer 
orbits  with  low  perigee,  high  apogee  and  eccentricity  near  0.7,  is  to  produce  an 
oscillation  in  perigee  height  with  an  amplitude  of  several  hundred  kilometres 
and  a period  of  several  years.  PROD  allows  the  forward  integration  of  the  orbit 
for  many  years,  and  the  date  when  the  perigee  is  driven  down  to  a height  near 
100  km  usually  indicates  the  likely  decay  date  to  within  a few  months.  In  quot- 
ing lifetimes,  we  usually  add  about  30  days  to  the  value  given  by  PROD,  to  allow 
for  a short  period  of  decay  under  air  drag,  which  is  often  found  to  occur  in 
practice.  For  satellites  of  long  life,  an  integration  interval  of  20  days  is 
often  used  with  PROD,  to  save  computer  time:  comparison  with  computations  having 
integration  intervals  of  10  days  and  5 days  usually  shows  no  great  differences. 
Shorter  intervals,  down  to  1 day,  are  used  when  more  detailed  results  are 

required.  Many  examples  of  orbital  evolution  calculated  by  means  of  PROD  have 
22-24 

been  given  , and  one  is  reprinted  as  Fig  19. 

There  are  two  main  difficulties  in  predicting  the  evolution  of  these  high- 
eccentricity  orbits  using  PROD.  The  first  problem  is  the  neglect  of  drag  in  the 
program,  and  this  is  discussed  in  section  5.3.  The  second  problem  is  the 
inevitable  inaccuracy  of  the  initial  values  of  perigee  height  used  in  the  computa- 
tion. The  USAF  NORAD  elements  are  nearly  always  the  only  elements  available  for 
these  satellites,  and  since  the  satellites  are  very  difficult  to  track,  often 
being  at  great  heights  in  the  northern  hemisphere,  the  orbits  are  often  based  on 
inadequate  numbers  of  observations  poorly  distributed  round  the  orbit,  so  that 
bias  errors  in  eccentricity  arise,  and  errors  of  up  to  30  km  or  even  50  km  in 
perigee  height  are  not  uncommon.  If  the  perigee  height  decreases  slowly  from 
200  km  to  100  km  at  the  end  of  the  life,  such  absolute  errors  in  initial  perigee 
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height,  which  tend  to  be  carried  through  the  integration,  can  obviously  cause 
serious  errors  in  the  estimates  of  decay  date. 

Despite  these  limitations,  PROD  is  often  successful  in  adverse  circum- 
stances. Fig  20  shows  a particularly  unfavourable  recent  example,  the  satellite 
1972-37A,  predicted  from  an  initial  orbit  in  March  1973.  The  initial  orbit  was 
probably  not  very  accurate,  and  the  integration  interval  was  rather  too  great 
for  good  accuracy  (25  days);  also  the  perigee  dipped  low  into  the  atmosphere  in 
November  1973,  and  decay  might  have  occurred  then,  but  did  not.  Yet  the  decay 
date  predicted  by  PROD  in  1973,  with  the  30-day  addition  mentioned  earlier,  was 
21  March  1977:  the  actual  decay  date  was  22  March  1977.  Fig  20  also  shows  the 
perigee  height,  calculated  from  the  decay  rate  n in  NORAD  bulletins,  as 
described  in  section  5.2,  for  the  occasions  when  n is  large  enough  to  be 
reliable. 


5.2  Calculation  of  perigee  height  from  decay  rate 

The  normal  method  of  calculating  perigee  height,  from  the  value  of 
a(l  - e) , is,  as  already  mentioned,  subject  to  errors  of  30  or  50  km  for  high- 
eccentricity  orbits.  If  the  perigee  height  is  near  200  km,  a change  of  30  km 
in  perigee  height  changes  the  air  density  at  perigee  and  hence  the  decay  rate 
by  a factor  of  about  3.  Such  errors  in  perigee  height  are  obviously  unacceptable 
if  accurate  forecasts  of  lifetime  are  to  be  made,  and  the  sensitivity  of  n to 
perigee  height  means  that  a better  estimate  of  perigee  height  can  often  be 
obtained  from  n if  the  mass/area  ratio  m/S  of  the  satellite  is  known.  If 
m/S  and  n are  both  known  accurate  to  about  ±10%,  an  accuracy  of  ±5  km  in 
perigee  height  can  be  achieved  for  heights  near  200  km,  or  an  accuracy  of  ±3  km 
if  perigee  height  is  near  150  km. 


The  relation  between  n and 


p may  be  found  from  equation  (22),  and 
P 


since  we  may  assume  that  the  perigee  height  is  200  ± 70  km,  equation  (23)  gives 
G(y^)  = 10  Pp*^H  , whence 
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10  p /h  1(e) 
P 


(26) 


Fig  21  gives  values  of  Z(e)  for  e 0.2,  and  Fig  22  gives  values  of  10  p^v^ 

from  CIRA  1972  for  100  < y < 260  km,  for  low  and  high  solar  activity  (exospheric 

P 

temperature  800  and  1100  K) . 

Fig  23  shows  the  variation  of  n(m/S)  with  y^  for  selected  values  of  e 
close  to  those  most  often  used:  for  example,  many  satellites  (including  the 
Molniyas)  have  eccentricities  between  0.7  and  0.75  during  nearly  the  whole 
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lifetime,  and  curves  are  given  for  these  two  values  but  not  for  e ■ 0.6,  a 
value  which  rarely  arises. 


To  use  Fig  23,  the  value  of  n (rev/day  ) is  calculated  from  the  two 

latest  USAF  bulletins  available  (or,  less  accurately,  from  the  value  of  n/2 

given  on  the  latest  bulletin),  and  multiplied  by  m/S.  The  corresponding  value 

of  y is  then  read  from  Fig  23,  for  the  appropriate  values  of  e and  solar 
P t t 

activity.  An  alternative  and  more  accurate  method  is  to  use  equation  (26) 

directly:  n(m/S)  is  calculated  as  before  and  divided  by  E(e)  , from  Fig  21, 

to  give  a value  of  lO^p^v^H  and  hence  y^  from  Fig  22. 


If  we  assume  that  the  value  of  y derived  from  n is  correct  and  call 

P 

it  y • , and  if  the  initial  value  of  y in  a PROD  run  derived  from  a(l  - e) 
pn  P 

is  y , then  PROD  should  be  re-run  with  e adjusted  to  give  the  correct 

perigee  height;  that  is,  e should  be  increased  by  Ay^/a  , where 

Ay  = y - y • . 
p pe  pn 


5.3  The  PTDEC  program 

The  program  PTDEC,  which  allows  for  drag  as  well  as  lunisolar  perturbations, 
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has  recently  been  developed  by  Palmer  , as  a shortened  version  of  the  program 
26 

POINT  (program  for  orbit  integration) . The  name  PTDEC  may  be  regarded  as  an 

acronym  for  POINT  to  decay.  The  integration  step  in  PTDEC  is  normally  1/96  of 

a revolution,  and  the  reference  atmosphere  used  is  a simplified  version  of  the 
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Jacchia  1965  model  , a precursor  of  CIRA  1972.  The  integration  ends  when  the 
perigee  height  drops  to  90  km  (or  any  specified  greater  value) . Values  of  S/m 
and  solar  10.7  cm  radiation  energy  have  to  be  supplied  by  the  user,  and  these 
values  are  kept  constant  during  the  computer  run. 

PTDEC  needs  much  more  computer  time  than  PROD,  and  is  essentially  a method 
for  use  in  the  last  few  months  of  a satellite's  life.  The  main  difficulty  is 
again  the  inaccuracy  in  the  initial  perigee  height,  and  it  is  usually  necessary 
to  have  a 'trial  run'  first,  with  e somewhat  reduced  to  avoid  the  possibility 
of  an  initial  perigee  height  lower  than  90  km,  and  then  a re-run  with  a corrected 
value  of  e , as  described  in  section  5.2. 

Fig  24  shows  the  variation  of  apogee  and  perigee  height  as  given  by  PTDEC 
for  the  last  two  months  in  the  life  of  1972-37A,  a typical  rapid-decaying 
satellite.  In  this  example  the  orbit  was  being  re-run  with  the  perigee  height 
given  by  PTDEC  arranged  to  coincide  with  the  perigee  height  deduced  from  the 
value  of  n in  NORAD  bulletin  329,  namely  207  km;  values  of  perigee  height 
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obtained  from  other  NORAD  bulletins  are  shown  as  circles.  The  lifetime  indi- 
cated by  PTDEC  was  63  days  from  January  25,  a date  at  which  the  semi-annaul 
correction  V/V  is  0.87,  from  Fig  17.  Thus  the  lifetime  given  by  PTDEC, 
corrected  for  the  semi-annual  variation,  is  55  days:  the  actual  lifetime  was 
56  days,  so  the  agreement  is  excellent.  In  running  PTDEC,  the  solar  F 

1 u • / 

index  was  taken  as  80,  which  proved  to  be  a good  approximation  to  the  actual 

value  in  February-March  1977,  namely  79.  The  value  of  S/m  was  taken  as 
2 

0.01  m /kg  in  PTDEC,  and  it  may  seem  surprising  that  such  an  accurate  lifetime 
estimate  was  obtained,  using  such  a 'round-number'  value  of  S/m  : but  since 
the  perigee  height  is  obtained  from  n assuming  this  same  value  of  S/m  , the 
initial  decay  rate  on  PTDEC  should  be  the  same  as  on  the  actual  orbit,  even  if 
S/m  is  in  error,  and  the  effects  of  errors  in  the  value  of  S/m  are  much  less 
than  might  be  expected. 

Fig  24  is  also  of  interest  in  its  own  right,  as  an  example  of  very  rapid 
decay.  Most  of  the  decrease  in  perigee  height  during  February  is  due  to  luni- 
solar  perturbations:  the  decrease  due  to  drag,  approximately  $H  lnCeg/e)  , is 
less  than  80  m (with  e^  = 0.750  and  e = 0.743  on  5 March).  The  decrease  in 
perigee  height  due  to  drag  only  becomes  appreciable  in  the  last  10  days  of  the 
life.  The  apogee  height  indicated  by  PTDEC  in  Fig  24  decreases  fairly  steadily 
during  the  week  before  decay  at  the  extremely  rapid  rate  of  3000  km  per  day. 

The  maximum  daily  decrease  is  3500  km,  two  days  before  decay. 

PTDEC  can,  of  course,  also  be  used  for  orbits  of  low  eccentricity,  and  it 

has  the  advantage  over  the  graphical  methods  that  it  automatically  takes  account 

of  the  variations  in  perigee  height  discussed  in  section  4.2.  However,  PTDEC 

does  not  allow  for  the  day-to-night  or  semi-annual  variations  in  density,  or  the 

future  variations  of  solar  activity;  so  these  still  have  to  be  estimated 

separately.  For  low-eccentricity  orbits  the  value  of  perigee  height  from  NORAD 
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elements,  with  the  appropriate  corrections  , should  be  accurate  to  2 km,  so 

there  is  no  need  to  obtain  perigee  height  from  n . Instead  it  is  possible  to 

use  n to  determine  S/m  for  use  in  the  program  if,  as  often  happens,  the 

value  of  S/m  is  not  known.  The  procedure  is  to  guess  S/m  for  a first  run, 

and  then  to  adjust  it  so  that  the  initial  value  of  n given  by  PTDEC  is  the 

same  as  that  of  the  actual  satellite.  The  values  of  m/S  for  a representative 

selection  of  actual  satellites  are  shown  in  Fig  15.  Values  of  m/S  near 
2 2 

100  kg/m  are  common,  and  0.01  m /kg  is  an  appropriate  value  to  use  as  a first 
guess  in  the  PTDEC  run. 


5.4  Decrease  in  drag  coefficient  at  low  altitudes 


The  theory  so  far  used  here  has  relied  on  the  assumption  that  the  drag 

coefficient  CQ  of  the  decaying  satellite  remains  constant,  and  when  a numerical 

value  is  needed,  as  in  equation  (20),  C has  been  taken  as  2.2.  This  is  the 
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value  appropriate  in  f ree-molecule  flow,  when  the  mean  free  path  of  the  mole- 
cules greatly  exceeds  the  dimensions  of  the  satellite.  When  the  mean  free  path 

becomes  comparable  to  the  dimensions  of  the  satellite,  however,  the  drag 
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coefficient  decreases  to  about  half  its  value  m f ree-molecule  flow  . The  mean 
free  path  (Fig  25)  decreases  from  10  ra  at  130  km  height  to  3 m at  120  km  height: 
so,  for  all  satellites  except  the  smallest,  the  remaining  lifetime  after  perigee 
has  descended  to  120  km  is  likely  to  be  double  that  expected  on  the  basis  of  a 
constant  value  of  . For  low-eccentricity  orbits,  this  effect  is  usually 
negligible  in  predicting  decay  dates  weeks  or  months  ahead,  because  perigee  does 
not  drop  below  120  km  until  within  a few  revolutions  of  decay.  The  effect  is 
also  negligible  for  high-eccentricity  orbits,  if  the  perigee  height  is  being 
driven  down  by  lunisolar  perturbations.  For  some  high-eccentricity  orbits,  how- 
ever, the  lunisolar  perturbations  can,  on  occasions  when  they  are  tending  to 
push  the  perigee  height  upwards,  combat  and  sometimes  outdo  the  decrease  in 
perigee  height  due  to  air  drag.  When  this  happens,  the  satellite  can  be  kept  in 

orbit  for  days  or  even  weeks  with  perigee  below  120  km.  An  example  is  1970-1 14F, 
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where  the  perigee  height  was  below  1 1 3 km  during  the  last  20  days  of  the  life  , 

and  the  drag  coefficient  was  calculated  to  be  0.85  ± 0.20. 

The  general  effect  of  the  decrease  in  when  the  mean  free  path 

decreases  to  a value  comparable  to  the  dimensions  of  the  satellite,  is  shown  in 

Fig  26,  which  gives  lifetimes  for  satellites  of  length  10  m (unbroken  lines)  and 

2 

1 m (broken  lines),  for  m/S  = 100  kg/m  . The  lifetime  is  proportional  to  m/S  , 
in  the  absence  of  strong  lunisolar  perturbations.  From  Fig  26  it  is  obvious 
that  the  change  in  has  a large  relative  effect;  but  most  satellites  have 

eccentricities  less  than  0.01  by  the  time  perigee  height  has  decreased  to  130  km, 
so  the  change  in  lifetime  for  these  satellites  is  less  than  0.2  day  (for 
m/S  = 100  kg/m^). 

For  the  few  satellites  which  are  affected  seriously  by  the  decrease  in 
CQ  , the  easiest  way  of  allowing  for  it  is  to  double  the  remaining  lifetime 
given  by  PTDEC  after  the  time,  t*  say,  when  the  perigee  first  descends  to  the 
height  where  the  mean  free  path  equals  the  length  of  the  satellite.  A more 
accurate  procedure,  necessary  when  lunisolar  perturbations  are  still  strong 
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(if  e > 0.5),  is  to  re-run  PTDEC  from  time  t*  onwards  with  the  original  value 
of  S/m  halved:  this  is  equivalent  to  reducing  to  half  its  value  in  free- 

molecule  flow. 

6 CONCLUDING  REMARKS 

The  aim  of  this  paper  has  been  to  outline  the  many  problems  encountered  in 
satellite  decay  predictions,  and  to  provide  solutions  to  as  many  as  possible  - 
but  regrettably  not  all. 

To  summarize,  the  lifetime  L of  satellites  in  orbits  of  eccentricities 
between  0 and  0.8  is  written  as  Q/n  , where  n is  the  current  decay  rate,  and 
values  of  Q are  given  first  in  Fig  2,  and  then  more  accurately  in  Figs  5,  6,  8 
and  10,  for  an  idealized  sperically  symmetrical  Earth  and  atmosphere.  Correc- 
tions to  these  values  for  the  effects  of  variations  in  perigee  height  and  semi- 
annual density  variations  are  shown  in  Figs  11-12  and  17-18  respectively,  and 
other  corrections  are  discussed  in  section  4.  For  orbits  which  suffer  serious 
lunisolar  perturbations,  numerical-integration  methods  are  needed,  and  the  use 
of  these  is  described  in  section  5. 

For  the  majority  of  non-manoeuvring  satellites,  these  methods  allow  the 
prediction  of  decay  dates  up  to  about  one  year  ahead  with  an  accuracy  of  about 
10%  of  the  remaining  lifetime,  if  sufficient  care  is  taken.  But  the  lifetime 
depends  on  a multiplicity  of  independent  variables,  including  the  rotational 
behaviour  of  the  satellite  (which  affects  the  cross-sectional  area),  the  time 
variations  of  the  atmosphere  (particularly  the  semi-annual  variation),  the 
future  variations  in  solar  activity  (for  which  prediction  is  poor  at  present), 
the  gravitational  fields  of  the  Earth,  Sun  and  Moon  (whose  effects  on  the  orbit 
themselves  depend  on  the  future  evolution  of  the  orbit) , and  the  future  values 
of  the  local  time  and  latitude  at  perigee  (which  depend  on  most  of  the  factors 
previously  mentioned).  Further  difficulties  are  that  current  orbital  elements 
may  not  be  available,  or  that  the  available  elements  are  inaccurate.  For  a 
small  minority  of  satellites  one  or  more  of  these  difficulties  arises  in  acute 
form,  and  can  greatly  degrade  the  accuracy  of  the  lifetime  estimates. 
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